Atomic three-body loss as a dynamical three-body interaction 



A. J. Daley/ J. M. Taylor,^ S. Diehl,i M. Baranov/ and P. ZoUeri 

'institute for Theoretical Physics, University of Innsbruck, A-6020 Innsbruck, Austria 
and Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, A-6020 Innsbruck, Austria 
'^Department of Physics, Massachusetts Institute of Technology, Building 6C-411, Cambridge, MA 02139, USA 

(Dated: 28 October 2008) 

We discuss how large three-body loss of atoms in an optical lattice can give rise to effective 
hard-core three-body interactions. For bosons, in addition to the usual atomic superfluid, a dimer 
superfluid can then be observed for attractive two-body interactions. The non-equilibrium dynamics 
of preparation and stability of these phases are studied in ID by combining time-dependent Density 
Matrix Renormalisation Group techniques with a quantum trajectories method. 



Cold atomic gases in optical lattices have proven a 
test bed for understanding novel quantum phases [Tj and 
non-equilibrium many-body dynamics [U Recently, 
Syassen et al. [4 showed that a strong two-body loss pro- 
cess for molecules in an optical lattice pLj could produce 
an effective, elastic hard-core repulsion and thus a Tonks 
gas [U This is related to the quantum Zeno effect: 
a large loss dynamically suppresses processes creating 
two-body occupation on a particular site. Whilst elastic 
two-body interactions occur in many systems, regimes 
where elastic 3-body interactions dominate are rare in 
nature. Here we discuss how the ubiquitous, though nor- 
mally undesirable three-body losses of atomic physics ex- 
periments can induce effective three-body interactions. 
These are associated with interesting quantum phases, 
including Pfafhan states [6J, and could be used to sta- 
bilise three-component Fermi mixtures [7\, assisting in 
the production of a colour superfluid state [8J . We inves- 
tigate Bosons in an optical lattice, where a three-body 
hard-core constraint stabilises the system with attrac- 
tive two-body interactions, and a dimer superfluid phase 
emerges. We focus on the dynamics of this intrinsically 
time-dependent system, both testing the hard-core con- 
straint for finite loss rates, and studying non-equilibrium 
properties including decay. In ID, the exact evolution is 
computed by combining time-dependent density matrix 
renormalisation group methods (t-DMRG) [3J [3] with a 
quantum trajectories approach from quantum optics [H]. 

Three-body recombination [10] in an optical lattice 
corresponds to decay into the continuum of unbound 
states, and thus loss from the lattice. This can be de- 
scribed by a master equation in the Markov approxima- 
tion [S] , which for atoms in the lowest band of an optical 
lattice can be projected onto the corresponding basis of 
Wannier functions [U |5j , associated with bosonic annihi- 
lation operators bi on site i. We can separate the master 
equation into terms which conserve particle number, cor- 
responding to an effective Hamiltonian H^s, and terms 
which remove three particles on a site: 
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FIG. 1: (a) Bosons in an optical lattice in the presence of 
three-body loss at a rate 73. (b) Example model parame- 
ters estimated for Gs at a magnetic field of 15 Gauss (where 
the scattering length is 20ao and the recombination length 
~ 500ao [lO] where ao is the Bohr radius) as a function of 
lattice depth Vb, showing 73 (solid line), U (dashed), and J 
(dotted). Values of 73 are obtained by integrating the mea- 
sured three-body recombination rates in free space over a 
state with three particles in a single Wannier function, (c) 
The probability that at least one loss event has occurred at 
time tJ = 2, beginning with a single particle on each of 10 
sites, and U/J = 3 (solid line), 5 (dashed), and 10 (dotted), 
computed using t-DMRG methods (see text for details). 

where p^"^ denotes the system density operator with n 
atoms and hi = blbi. The dominant loss term is on-site 
three-body decay [H] and 73 is the corresponding rate. 
The effective Hamiltonian is 

- H--^^T.(birbi (1) 

i 

H = -J^b\bj + ^^h,{h,-l) + ^eih„{2) 

with J the nearest neighbor tunneling amplitude, U the 
elastic two-body interaction, and et the local potential. 
The Hamiltonian is valid in the limit where J, £i, Un, <C 
u, with LO the band gap and n the mean density. In 
an experiment these parameters, in particular U, may 
be tuned whilst 73 remains constant and large [15]. In 
Fig. Ic we show example values of 73, U, and J using 
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numbers for Caesium as a function of lattice depth. 

If we begin in a pure state with N particles, then loss 
processes lead to heating, in that they produce a mixed 
state of different particle numbers. Within a fixed par- 
ticle number sector, the dynamics are described by -ffcff • 
Three-body interactions emerge most clearly in the limit 
of rapid decay: 73 J, U, e^. If we define the projector P 
onto the subspace of states with at most two atoms per 
site and Q = 1 — P, then in second order perturbation 
theory we obtain the effective model 

2i 6 
HZ. w PHP PHQHP = PHP -I PS^ c\cjP, 

73 ^3 i 

(3) 

where Cj = {b']/V^)Y. fegjv ^k, and Nj denotes the set of 
nearest neighbours of site j. The term PHP describes 
the Hubbard dynamics, Eq. ([2]), supplemented by the 
hard-core constraint {b\Y — 0. Furthermore, the effec- 
tive loss rates decrease as J^/73 [H]. 

Thus, we see the clear emergence of a three-body hard- 
core constraint in the limit 73/ J ^ 1. We can study 
the physics of the projected model PHP to obtain a 
qualitative understanding of the quantum phases asso- 
ciated with the projection. However, the residual loss 
processes make this system intrinsically time-dependent, 
and can give rise to heating. We therefore study the full 
non-equilibrium dynamics, by combining t-DMRG meth- 
ods [21 [3] with an expression of the master equation as 
an average over quantum trajectories [3]. Each stochas- 
tic trajectory begins from an initial pure state (sampled 
from the initial density matrix), and can be interpreted 
as describing a single experimental run, in which losses 
occurred at particular times i„ and on sites i„. The evo- 
lution is described by the non-Hcrmitian i/cffi except for 
times tn, where losses (or quantum jumps) occur, 

in^^\m)-H..m)); i^(^>t)) = |g;j^;;j;|| , (4) 

where the jump operator Ci — &f corresponds to three- 
body loss on site i. In stochastic simulation of the master 
equation, the times tn are points where the norm of the 
state falls below a randomly chosen threshold. At these 
times, a random jump operator is selected according to 
the probabilities pi^ cx {ip{tn)\ClCi\4'{tn)) and applied. 
In this way we can both investigate individual trajec- 
tories and compute expectation values from the master 
equation. The latter is performed by stochastic average 
over both initial states and over jump events, which con- 
verges rapidly as the number of trajectories is increased. 

The need to simulate many trajectories for conver- 
gence is offset by the efficiency of simulating states rather 
than density matrices. In contrast to directly applying 
t-DMRG to density matrices [11], we can simulate the 
master equation without squaring the local Hilbert space 
dimension, and by time-evolving states with fixed parti- 
cle number, we also make use of existing optimisations 



for conserved quantities. Despite the application of lo- 
cal jump operators, we find the evolution quite efficient, 
especially for small numbers of jumps [17] . 

As an example of the suppression of loss, we consider 
preparing a homogenous initial state at unit filling in a 
deep optical lattice where U/J ^ 00. At time i = we 
suddenly ramp the lattice to a finite depth, and observe 
the probability p that a single three-body loss event has 
occurred as a function of time. In Fig. lb we plot this 
probability for different U/J as a function of 73/ J. We 
see a clear suppression of loss rates for large 73 / J, and 
also a substantial decrease for larger U / J , resulting from 
the decreased amplitude for doubly occupied sites. 

In the limit of large 73, it is instructive to study the 
equilibrium phase diagram of the projected Hamiltonian 
PHP. For U/J > 0, we observe the well-known Mott 
Insulator (MI) and atomic superfiuid phases of the Bose- 
Hubbard model. However, the three-body hard-core con- 
dition will also stabilise the system for U/J < 0, where 
we find a dimer superfiuid phase (see Fig. 2a). This is 
characterised by the vanishing of the order parameter 
signalling superfluidity of single atoms (ASF) {{bi) = 0), 
while a dimer superfluidity (DSF) order parameter per- 
sists ((6f) ^ 0). The superfluid regimes are connected via 
a quantum phase transition associated with the sponta- 
neous breaking of a discrete Z2 symmetry, reminiscent of 
an Ising transition [T^]: the DSF order parameter trans- 
forms with the double phase ~ exp2i6' compared to the 
ASF order parameter ~ expi^. Consequently, the sym- 
metry 6 ^ 6 + IT exhibited by the DSF order parameter 
is broken when reaching the ASF phase. 

We can obtain a qualitative mean-fleld picture using a 
homogeneous Gutzwiller ansatz wavefunction, given for 
the projected Hilbert space by I^P) = Yli where 
1^-)^ = roe*'^«|0)j -I- rie'-^ill)* + r2e''*2|2),. Normalisation 
implies J^a = 1; while the filling is n = rf + 2r| < 2. 
To examine the phases, we find the energy E/M'^ = 
(^-liJl^-), E{ra,(j)a) = Urj - Jzrl[rl + 2V2r2roCOS^ + 
2r|], where ^ — <j>2 + 4>() ^ 2(/)i and M'' is the number 
of lattice sites. For any the energy is minimized for 
<i> an integer multiple of 27r. For ri,ro ^ (placing us 
in the ASF phase with |(6,)p = r?(ro + V2r2)^ + 0), 
the phase-locking expression contributes a "source" term 
linear in r2 to the energy, and consequently the mini- 
mum of the energy cannot be located at r2 = 0. Thus, a 
finite atomic condensate always implies a dimer compo- 
nent = 2(ror2)^. However, the reverse is not true 
and one may have pure dimer superfluidity without an 
atomic condensate. 

At fixed n, the energy is a function, e.g., of r\ alone. 
For the second order transition found within our mean 
field theory, an instability for atomic superfluidity is 
indicated by its mass term crossing zero, d^E(r\ — 
0) /dr\ = 0. This leads to a critical interaction strength 
for the ASF-D SF transition, Uc/{Jz) = -2{l + n/2 + 
2yjn{l — n/2)). Within the DSF phase the order pa- 
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rameter obeys = ^(1 — n/2) independent of the 

interaction strength. For n — s- 2, we approach a MI state 
in a second order transition. At n = 1 we find that 
the ASF-DSF transition takes place at the same cou- 
phng strength as the ASF-MI transition, but with the 
opposite sign. The complete mean field phase diagram 
in the plane of density and interaction strength is plotted 
in Fig. 2a. 

From the last term in Eq. we can estimate the 
initial loss rate from the ground state Gutzwiller wave 
function. We obtain the rate jeS = 3 J^2:/73M'^((n^) — 
n){n + which is zero in the MI, and oc for 

the DSF, 7cff = 3J'^ zM'^n'^ /-f^. In the DSF phase, the 
critical temperature Tc oc ri^^^ at low densities, and 
the energy density deposited by a single loss, Ai?ioss = 
{z + l)\U\n/{2M'^). The number of independent loss 
events needed to melt the DSF is then proportional to 
Tc/AEiassi E^nd the melting time strongly decreases for in- 
creasing density, proportional to ^3/{\U\{z + Vjj'^zn'/^). 

These qualitative features are reproduced in ID, as 
supported by numerical calculation of the ground state 
for PHP. In Figs. 2b, c we show the characterisation of 
the crossover between the ASF and DSF regimes in ID 
via the off-diagonal elements of the single particle density 
matrix, S{i,j) = {b\bj), and the dimer density matrix 
D{i,j) = (blblbjbj). In the MI regime, the off-diagonal 
elements of S{i,j) and D(i,j) decay exponentially. As 
we enter the superfluid regime, quasi-long range order 
is visible in the polynomial decay (linear on the loga- 
rithmic scale). As U/J is made more negative, we see 
a return to exponential decay for the off-diagonal ele- 
ments of S{i,j), but the off-diagonal elements of D{i,j) 
still decay polynomially and, indeed, increase in mag- 
nitude. This characterises the DSF regime in ID. Here, 
the transition to the DSF and MI regimes occurs at much 
smaller \U/J\ than in higher dimensions, but these two 
transitions again occur at similar \U/J\ for n— I. 

A dimer superfluid phase could be prepared in an ex- 
periment beginning in states that have very small am- 
plitude of three-body occupation due to two-body inter- 
actions. Hamiltonian parameters could then be ramped 
adiabatically across the phase boundary, making use of 
dynamical suppression of three-body occupation. We 
study two such scenarios as illustrated in Fig. 3a: (i) 
We begin from a Mott Insulator state, and then ramp 
from J7/J = 30 to U/J = —8 in order to enter the dimer 
superfluid regime (which is intuitive, but associated with 
large probabihty of decay) ; or (ii) We apply a super lattice 
that raises the energy of certain sites and form a Mott 
Insulator with two particles per site in the lowest wells, 
then switch the interaction rapidly to U/J = —8 on a 
timescale much faster than tunnelling between the lowest 
wells and ramp down the superlattice. In each case, we 
compute dynamics using t-DMRG methods with quan- 
tum trajectories, and consider an on-site three-body loss 
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FIG. 2; Equilibrium analysis of the projected Bose-Hubbard 
model PHP. (a) Mean-field phase diagram as a func- 
tion of U/ (Jz) and density, n. (b,c) Magnitude of off- 
diagonal elements of (b) the single particle density ma- 
trix 15(1,^)1 = \{blbj)\ and (c) the dimer density matrix, 
\D{i,j)\ = \{blblbjbj)\, as a function of \i - j\, for U/J =10 
(thin solid line), 5 (dotted), (dot-dash), -5 (dashed), and -10 
(thick solid line). These results are computed for 20 particles 
on 20 lattice sites in ID, with box boundary conditions and 
i+ j — 21, using imaginary time-evolution in t-DMRG, and 
plotted on a logarithmic scale. 



with 73/J = 250. We ramp parameters on a timescale 
where without loss the ground state will be reached with 
minimal presence of excited states. 

In Figs. 3b, c we show the time dependence of the sum 
of kinetic and interaction energies, and of the total parti- 
cle number for example trajectories. For each ramp type, 
we choose a "lossless" trajectory, where the ground state 
is reached adiabatically, and a "lossy" trajectory, where 
three-body loss events lead to heating of the system (as 
holes are produced that correspond to excited states). 

In Fig. 3d, we compare the probability for each type of 
ramp of producing the "lossless" trajectory. We see that 
for the ramp from the MI state, where we pass through 
a region of small U/J, the probability of such a trajec- 
tory is essentially zero. For the superlattice ramp, on 
the other hand, it is much more likely that we will ob- 
tain the ground state from a randomly chosen trajectory. 
This is both because the superlattice allows us to use 
large \U/ J\, so that the two-body interaction reduces the 
probability of triply occupied sites, and because the su- 
perlattice facilitates the choice of a lower density in the 
system, which signiflcantly reduces the effective loss rate. 

In Fig. 4 we show the local density as a function of 
time and the final dimer density matrix D(i,j) for (a) 
the "lossless" and (b) the "lossy" trajectories of Fig. 3c. 
When a loss occurs it affects not just the density for the 
site on which it occurs, but also on neighbouring sites due 
to the knowledge that we obtain of the position of the 
remaining particles. We also see clearly the destruction 
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FIG. 3: Dynamics of adiabatic ramps into a dimer superfluid 
regime, (a) We begin with (i) a Mott-insulator state (ramping 
U / J), and (ii) a state with pre-prepared dimers in a superlat- 
tice (removing the superlattice). (b)-(c) The sum of kinetic 
{Ek) and interaction {Ei) energy and (inset) particle num- 
ber as a function of time for two example trajectories, one 
with no loss events (dashed lines) and one with several loss 
events (solid lines). Here, (b) shows a ramp from f//J = 30 
to U/J = -8, with U{t) = aJ/(100 + it J) + 7, with a and 7 
ramp parameters, and (c) shows a ramp with a superlattice 
potential, ei — Vo cos(27r//3), where Vo ~ 30Jexp(— O.lt J), 
adjusted so that Vo{tJ = 100) = 0, with fixed U/J = -8. In 
each case, 73 = 250 J. For (b), we use 20 atoms on 20 lattice 
sites, for (c), 14 atoms on 23 lattice sites, (d) Plot showing 
the probability that no loss event has occurred after time t 
for the ramps in (b) (dashed line) and (c) (solid). 



of correlations in the region of the system where the loss 
occurs. Note, however, that off-diagonal correlations are 
still visible in parts of the system, and that a single loss 
event does not always destroy the properties of the final 
state entirely. As discussed above, the probability of an 
individual loss event increases with system size, but a 
single loss event will change the character of the final 
state less. 

The three-body interactions discussed here could have 
applications to producing Pfaffian-like states, and also 
for fermions, especially stabilising three-component mix- 
tures. The combination of quantum trajectories methods 
with t-DMRG has potential applications in the simula- 
tion of other classes of master equations. In addition, 
there are open questions regarding the nature of the ASF- 
DSF phase transition, which we will address in more de- 
tail within a quantum field theoretical treatment of the 
attractive Bose-Hubbard model with a three-body hard- 
core constraint |13j . 
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FIG. 4: Comparison of (a) lossless and (b) lossy trajectories 
from Fig. 3c. We show the mean density (rii) as a function 
of position and time, and magnitude of the dimer correlation 
function \D{i,j)\ {i 7^ j) at the end of the ramp. 
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